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Analysis  of  Scalar  Datasets  on 
Multi-Resolution  Geometric  Models 

Alexandre  Gerussi  and  Georges-Pierre  Bonneau 


Abstract.  Recently,  multi-resolution  methods  based  on  non-nested 
spaces  were  introduced  to  allow  the  visualization  and  approximation  of 
functions  defined  on  irregular  triangulations  [3,4,5].  This  paper  comes  back 
to  these  methods  and  shows  more  precisely  how  the  subdivision/prediction 
/correction  scheme  of  ordinary  wavelet-based  multi-resolution  analysis 
(MRA)is  also  present  in  that  framework.  As  an  illustration,  it  is  demon- 
strated how  it  can  be  applied  in  two  of  the  classical  issues  of  MRA:  com- 
pression and  level-of-detail  editing.  We  also  show  that  the  framework 
can  be  used  for  the  analysis  and  approximation  of  scalar  data  defined  on 
meshes  with  arbitrary  topology,  thus  extending  our  previous  results  in  the 
plane  and  the  sphere.  Here  again,  the  link  with  the  corresponding  classical 
multi-resolution  scheme  of  [6]  as  well  as  decimation  methods  is  made. 


§1.  Introduction 

In  the  last  few  years,  the  problem  of  simplifying  huge  3D  triangular  meshes, 
for  the  purpose  of  e.g .,  visualization,  transmission  or  storage,  has  received 
considerable  attention.  Among  those  works,  two  major  approaches  can  be 
found.  In  the  case  of  regular  meshes,  the  use  of  a wavelet-based  framework 
has  proven  to  be  a powerful  solution  [6,9,14].  On  the  other  hand,  when  meshes 
are  not  regular,  the  approach  has  been  to  simplify  the  mesh  by  applying 
a sequence  of  elementary  geometric  simplification  operations,  such  as  vertex 
removals,  edge  collapses  or  triangle  collapses,  the  order  of  removal  being  driven 
by  a greedy  algorithm  [1,8,11],  We  refer  to  this  latter  approach  as  a decimation 
approach. 

This  paper  is  concerned  with  the  simplification  of  data  that  is  defined 
on  a surface  by  means  of  a triangulation.  This  topic  is  closely  related  to 
surface  simplification,  since  a triangular  mesh  can  be  seen,  at  least  locally,  as 
the  graph  of  a piecewise  linear  function  supported  by  a triangulation.  Here 
again,  when  the  surface  is  well-known,  for  example  a plane  square,  a sphere, 
a cylinder,  etc.,  several  wavelet  based  approaches  have  been  employed  [9,14], 
with  regular  underlying  meshes. 
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In  [3,4,5]  the  concept  of  non-nested  MRA  was  introduced  and  applied 
to  the  approximation  and  progressive  visualization  of  piecewise  constant  or 
linear  functions  defined  on  arbitrary  planar  or  spherical  meshes.  In  Section  2, 
we  investigate  the  relationship  between  the  non-nested  framework  and  irreg- 
ular subdivision.  Examples  of  compression  and  level-of-detail  editing  in  that 
framework  are  given  in  Section  3. 

Section  4 focuses  on  the  approximation  and  visualization  of  functions 
defined  on  triangular  meshes  with  arbitrary  topology.  Here,  an  additional 
difficulty  is  that  the  surface  which  supports  the  function  is  also  altered  by  the 
approximation  process.  Like  in  the  case  of  surface  simplification,  a wavelet 
approach  can  successfully  be  applied  when  the  original  mesh  has  subdivision 
connectivity  [6,12],  Otherwise,  in  an  irregular  setting,  a decimation  approach 
is  usually  employed,  and  the  function  is  approximated  during  the  simplifica- 
tion process  [1,11].  We  show  how  to  apply  our  framework  to  functions  defined 
on  such  general  meshes.  We  will  see  that  it  makes  the  link  between  the 
wavelet-based  approach  available  for  subdivision  surfaces  and  the  decimation 
model.  As  in  our  previous  papers,  the  function  in  its  multi-resolution  form 
is  described  by  a coarse  approximation  defined  on  the  simplified  mesh,  and 
a sequence  of  detail  coefficients  that  are  used  for  the  reconstruction  of  the 
function  on  every  LOD  up  to  the  original  mesh.  Our  scheme  is  fully  bijective: 
The  function  multi-resolution  representation  has  the  same  size  as  the  original 
one.  The  approximation  process  performs  Z-2  approximation  of  the  data,  but 
other  types  of  approximation  are  also  possible. 

§2.  Non-Nested  Framework  and  Irregular  Subdivision 

In  Section  2.1  we  briefly  review  the  non-nested  decomposition  scheme  de- 
scribed in  [3,5].  Section  2.2  makes  the  link  with  the  notion  of  irregular  subdi- 
vision. 

2.1.  Decomposition  scheme 

For  simplicity,  every  space  is  supposed  to  have  finite  dimension.  Let  ft  be  a 
measurable  domain  and  V',  i = Q,...,N,  a sequence  of  subspaces  of 
These  spaces  do  not  have  to  be  nested  but  will  in  general  be  “growing”  in  the 
sense  that  dim(W)  < dim(W+1).  Now  let  / = /n  be  a function  in  the  finest 
subspace  VN . In  classical  MRA,  the  spaces  are  nested  and  the  link  between 
V and  VN  is  made  by  taking  a complementary  space  WN~l  of  VN~l  in 
VN , that  is 

yN  = yN-l^wN-\ 

Now  if  we  write  /jv  = }n-\  + 9n-i  according  to  the  space  decomposition, 
/jv_ i can  be  seen  as  an  approximation  of  / ^ in  VN_1,  and  f/jv-i  as  the  detail 
needed  to  recover  the  original  function  from  its  approximation.  By  repeating 
the  decomposition,  one  obtains 

In  = fo  + 9o  H + 9n-i 

which  corresponds  to  the  space  decomposition  VN  = V°  © 1 W‘. 
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Notice  that  in  this  case  /,  = Py<  (/i+i),  where  Pyi  is  the  projector  on  V'  with 
direction  W\  We  return  to  the  general  case  and  suppose  that  a linear  “pro- 
jector” P1  : Vi+1  -»  V1  is  given.  To  avoid  technical  details,  these  projectors 
are  required  to  be  surjective,  but  the  results  in  this  paragraph  also  hold  if 
this  is  not  the  case.  Let  W*  be  the  kernel  of  P*,  and  V'1  be  a complementary 
space  of  W*  in  V1+1.  We  now  observe  that  the  restriction  of  Pl  to  V1  is  a 
bijective  operator,  having  the  same  range  as  P%.  Thus  if  fi  = P'(/j+i)  and 
9i  = Q'(/i+i),  where  Ql  is  the  projector  on  Wl  (defined  by  the  choice  of  V1), 
the  following  reconstruction  formula  holds: 


fi+i  = Inv(P|pi)(/«)  + 9i-  (1) 

Again,  by  iterating  this  decomposition,  we  obtain  a coarse  approximation 
fo  = P°  o • • • o Pjv~1(/at)  and  “detail”  functions  go, , Qn-i- 

We  now  take  a look  at  the  reconstruction  process.  Denote  by  Sl  the 
inverse  of  P.G,.,  Sl  : V1  — > V1  C Vj+i.  The  complete  coarse-to-fine  reconstruc- 
tion formula  is  obtained  by  iterating  the  reconstruction  formula  (1): 

}N  = sN~1o- ■■oS°{f0)  + SN~1o- ■■oS1(g0)  + --  ■ + SN~1(gN-2)  + 9N-i-  (2) 


2.2.  Approximating  spaces,  subdivision  spaces  and  scaling  spaces 

We  are  going  to  see  that  the  previous  scheme  can  actually  be  considered  in 
two  different  ways.  Until  now,  it  was  implicitly  assumed  that  the  spaces 
V1  were  playing  the  role  of  the  scaling  spaces  in  classical  MRA.  Under  this 
assumption,  we  conceptually  have  a really  non-nested  framework;  if  the  spaces 
were  nested,  the  operators  Sl  would  be  the  identity  (formally  injecting  V1  in 
Vt+1).  However,  we  will  not  call  them  scaling  spaces  but  approximation  spaces, 
and  keep  the  term  “scaling”  for  other  spaces  that  are  going  to  be  defined  below. 

In  the  non-nested  framework,  one  loses  the  notion  of  subdivision  (or  cas- 
cade algorithm).  However,  looking  at  things  slightly  differently  allows  subdi- 
vision to  fit  into  the  non-nested  scheme.  To  show  this,  consider  the  operators 
S’  as  subdivision  operators,  and  call  the  spaces  V'  subdivision  spaces  accord- 
ingly. This  means  that  we  start  from  a function  fo  E V°  and  iteratively 
subdivide  it  into  fi,  fi,-  ■ • fiv  using  the  formula  fi+i  = Sl(fi).  The  notion 
of  subdivision  used  here  is  very  general.  In  that  context,  classical  regular 
or  semi-regular  subdivision  schemes  would  give  rise  to  a nested  sequence  of 
subdivision  spaces.  But  completely  irregular  schemes  would  require  the  non- 
nested framework  to  be  fitted  in.  The  use  of  non-nested  MRA  was  introduced 
in  [2],  and  was  later  applied  to  triangular  schemes  in  [3,4,5].  Recently,  another 
approach  on  general  irregular  schemes  was  proposed  in  the  work  of  Sweldens 
and  Guskov  [10,13]. 

We  now  define  the  scaling  spaces.  Like  in  classical  MRA,  they  are  the 
spaces  containing  the  limit  functions  resulting  from  the  subdivision  process. 
Since  very  little  is  known  about  the  convergence  of  such  schemes,  we  won’t 
push  the  subdivision  to  infinity,  but  restrict  ourselves  to  an  integer  N.  This 
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makes  sense  since,  when  using  such  a framework,  one  usually  starts  from  an 
initial  triangulation  Tjv  which  is  coarsened  to  Tq.  This  contrasts  with  the 
traditional  approach  in  subdivision  where  one  starts  from  a base  mesh  and 
subdivides  it  according  to  a systematic  rule. 

The  scaling  spaces,  for  i — 0, . . . , N,  are  defined  as 

V''N  = {/■'"  | fk  € V% 

where  f'k’N  = S1,N(fk),  and  S',N  = SN_1  o ■■■  o S'.  This  operator  carries 
functions  in  V'  through  N — i subdivision  steps  to  functions  in  V',N . 

Notice  that,  for  every  i,  V',N  C VN  and  moreover,  that  V',N  C V'+1,N . 
Now  fix  a basis  (</>*)  for  V'  and  (ip')  for  W'.  If  (a*)  and  (bl)  denote  the  coor- 
dinates of  fi  € V'  and  gr  £ W'  with  respect  to  these  basis,  the  reconstruction 
formula  can  be  re-written  as 


}n  = Y,aitiN  + Y Y^Y* 

k i—j  k 

for  any  j = 0, . . . , N,  which  corresponds  to  the  decomposition  of  //v  in 
VN,N  = VN  (using  the  same  notations  as  above:  ip*kN  = Si'N(tpk),  and 

i’Y  = si+1'N(iPi)). 

Although  formally  identical,  considering  the  V’s  as  approximation  or 
subdivision  spaces  changes  the  aspect  of  several  questions.  For  example,  in 
the  problem  of  error  measure  in  the  context  of  approximation  spaces,  we 
are  interested  in  ||//v  — fi  ||,  whereas  in  the  other  context  we  are  looking  for 

UN-fYi 

§3.  Application  to  Data  Compression  and  LOD  Editing 

In  this  section  we  show  how  the  framework  can  be  used  on  functions  defined 
over  irregular  triangulations  to  achieve  data  compression  and  LOD  editing, 
which  are  both  standard  applications  of  MRA.  The  context  here  is  the  planar 
or  spherical  setting  of  [4,5],  only  the  filters  need  to  be  changed.  Indeed,  the 
analysis  operator  used  in  those  papers  was  the  orthogonal  projector  (the  goal 
being  progressive  visualization).  However,  for  compression  it  is  often  useful  to 
know  in  advance  the  error  between  the  original  function  and  its  approximation, 
in  terms  of  the  wavelet  coefficients  that  were  used  in  the  reconstruction.  This 
will  be  achieved  by  designing  new  filters. 

3.1.  Isometric  subdivision 

In  order  for  the  detail  coefficients  to  have  the  error  measure  property,  the 
synthesis  operator  is  required  to  be  an  isometry.  Indeed,  suppose  that 

• S'  : V'  — > V'+1  is  an  isometric  operator,  Vi  = 0, . . . , N — 1, 

• the  complementary  spaces  V'*  are  chosen  orthogonal  to  the  W' , Vi  = 
0,  • • • , AT  — 1. 
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When  the  latter  condition  is  fulfilled,  we  say  that  we  are  in  a semi-orthogonal 
framework.  If  /jv  denotes  the  original  function,  then,  according  to  the  global 
reconstruction  formula  (2),  the  quantity  ||/v  — /q’"||2  which  measures  the 
contribution  of  the  correction  steps  in  the  reconstruction  process  is 

IIS"-1  o • • • O S\go)  + ■ • ■ + SN~1(gN-i)  + ^-i||2- 

Because  of  semi-orthogonality,  it  is  equal  to 

IIS"-1  o • • • o S1(g0)  + ■■■  + S"-W2)||2  + ll^-ill2, 

and  because  S"-1  is  an  isometric  operator,  we  can  factorize  and  remove  it 
from  the  first  term  above,  and  then  iterate  the  operation  to  get 


ll/iv-/o0’"ll2=£W 

1=0 

Notice  that  even  if  the  subdivision  operators  are  not  isometric, 

II fN  - f°o’N\\2  < IIS"-1!!2  ■ ■ • ||S1||2||5o||2  + • ■ ■ + ||S"-1||2||3n_2||2  + ||ffiv-i||2 

still  holds  in  the  semi-orthogonal  setting.  Let  fN  denote  the  partially  re- 
constructed function,  and  suppose  in  addition  that  we  are  in  an  orthonormal 
framework,  that  is,  the  functions  if'  form  an  orthonormal  basis  of  W'.  Pro- 
ceeding as  above  leads  to 


ii/n-/*ii2=  ££404)2, 

i=0  k 

where  the  6J.’s  are  the  wavelet  coefficients  and  e'k  equals  1 whenever  b\  is  taken 
in  the  reconstruction  and  0 otherwise.  Consequently,  in  this  setting  we  have 
an  error  measure  in  terms  of  wavelet  coefficients. 

3.2.  An  isometric  subdivision  operator 

The  idea  behind  this  construction  is  that  the  corresponding  analysis  oper- 
ator should  have  reasonable  approximation  quality,  which  seems  intuitively 
required  to  achieve  compression.  Accordingly,  the  projection  operator  used 
in  [4,5]  is  taken  as  a starting  point,  the  problem  being  to  approximate  it  by 
means  of  an  isometry.  Let  P'  \ P,+1  — > V*  be  that  operator.  The  first  step 
is  to  find  the  matrix  of  Pl  with  respect  to  some  orthonormal  basis  (eI+1)  and 
(el)  of  Vl+1  and  V1.  Now  let  UDV  be  the  singular  value  decomposition  of 
that  matrix.  This  decomposition  admits  the  following  interpretation: 

• V is  the  matrix  of  an  isometric  operator  of  P‘+1  in  the  basis  (eI+1)  since 
it  is  unitary  and  the  basis  is  orthonormal. 
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• D is  a £ x c matrix  (£  < c)  whose  diagonal  coefficients  are  all  non-negative 
and  < 1.  The  diagonal  matrix  formed  by  the  positive  entries  of  D is  the 
matrix  of  a bijective  operator  mapping  V1  onto  Ran(_P‘). 

• Like  V , the  matrix  U is  an  isometric  operator  of  V 1 . 

Let  D be  the  matrix  obtained  by  replacing  every  positive  diagonal  element  of 
D by  1.  This  amounts  to  turning  the  bijective  operator  above  into  an  isometric 
one,  and  thus  JJ DV  is  the  matrix  of  an  operator  Pl  whose  corresponding 
subdivision  operator  S'  is  isometric.  A few  remarks  can  be  made  about  this 
construction: 

• P'  does  not  depend  on  a particular  choice  of  orthonormal  basis. 

• Pl  is  not  the  best  isometric  approximation  of  P*  in  terms  of  L 2 norm 
of  operators,  but  it  can  be  shown  that  it  is  the  best  with  respect  to  the 
Frobenius  norm. 

• The  diagonal  coefficient  of  D are  by  definition  the  cosines  of  the  angles 
between  the  spaces  V 1 and  V'+1  (see,  e.g.,  [7]  Chapter  1).  In  the  nested 
case,  they  would  all  be  equal  to  1,  and  the  corresponding  subdivision 
operator  would  be  the  identity. 

• This  method  could  be  used  to  approximate  operators  by  means  of  simi- 
larities, by  replacing  the  entries  of  D by  an  appropriate  scalar  a instead 
of  1.  Although  better  in  terms  of  approximation  quality,  this  leads  to  bad 
visual  results  since  the  resulting  subdivision  operator  doesn’t  reproduce 
constants  if  a ^ 1. 

Notice  that  in  the  context  of  [4,5],  this  approximation  is  always  computed 
locally,  leading  to  a global  algorithm  in  linear  time. 

3.3.  Examples 

As  it  is  mentioned  in  the  beginning  of  this  section,  these  examples  were  created 
using  the  setting  described  in  [3,4,5] ; the  reader  is  invited  to  look  there  for 
details.  The  initial  triangulation  is  completely  irregular,  generated  by  random 
vertex  insertion.  In  Figure  1,  a piecewise  linear  setting  is  used  for  LOD  editing. 
The  function  is  edited  at  a coarse  resolution,  by  pulling  values  up  (— > white) 
at  some  vertices,  and  then  adding  detail  coefficients  back.  Figure  2 shows  an 
example  of  data  compression  when  the  approximation  spaces  are  spaces  of 
piecewise  constant  functions  over  triangulations  generated  from  the  original 
one  by  homogeneous  decimation.  This  last  setting  is  the  full  generalization  of 
Haar  wavelets  to  irregular  grids,  as  it  would  lead  to  them  in  the  regular  case. 

§4.  Scalar  Datasets  on  Irregular  Meshes  with  Arbitrary  Topology 

In  this  section  we  describe  how  the  non-nested  framework  can  be  used  to 
handle  scalar  attributes  defined  on  meshes  with  any  topology.  As  a starting 
point,  we  take  a multi-resolution  decimation  model,  based  on  the  vertex- 
removal  (VR)  operation  to  simplify  the  geometry.  This  means  we  assume 
that  an  initial  fine  mesh  is  given  along  with  its  associated  sequence  of  VR’s 
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iii.  50  points  edited  iv.  Reconstruction 


Fig.  1.  Picture  design  using  level-of-detail  editing. 


Original:  50000  triangles  Compressed  to  5%,  6%  error 

Fig.  2.  Compression  of  a piece-wise  constant  function. 

that  can  be  progressively  applied  to  decimate  the  mesh  down  to  a base  mesh. 
In  addition,  we  suppose  that  the  scalar  attributes  are  defined  by  means  of 
piecewise  constant  and/or  linear  functions  parameterized  on  the  initial  mesh. 
In  what  follows,  PI(V)  denotes  the  polygon  of  influence  of  a vertex  V;  it  is  the 
polygonal  area  delimited  by  the  1-neighbours  of  V.  In  order  to  apply  the  non- 
nested framework,  approximation  spaces  and  approximation  operators  need 
to  be  defined. 

4.1.  Local  mapping 

Let  M',i  = No,  denote  the  triangular  mesh  consisting  of  i vertices  (the 

original  mesh  after  N — i VR  operations).  Let  T(Ml)  be  the  space  of  real- 
valued functions  defined  on  M\  and  Cl  (resp.  £*)  be  the  subspace  of  functions 
oiF{Ml)  that  are  piecewise  constant  (resp.  linear)  on  each  triangle  of  M*.  We 
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refer  to  these  spaces  as  the  upper  approximation  spaces.  Each  VR  alters  locally 
the  surface,  thus  functions  of  T(Ml+l ) and  T{Ml)  are  defined  on  different 
domains.  To  define  an  approximation  problem,  a common  parameterization 
for  these  functions  is  required.  To  that  end,  we  assume  that  for  every  VR  of 
vertex  Vi,  a local  one  to  one  projection  II1  of  PI(Vt)  onto  a plane  is  also  known. 
The  reader  can  refer  to  [8]  for  a study  of  the  existence  and  determination  of 
such  a projection.  We  use  II'  to  consider  the  change  of  parameterization 
Hl  : M'  — > M'+1  as  the  mapping  defined  by  IP  o H'  = IP  over  PI(Vi)  and  by 
the  identity  outside. 

4.2.  Scaling  spaces  and  data  decomposition 

Let  K.  stand  for  C or  C.  To  a function  / € K.1  we  associate  the  function 
f i,i-i  _ an(j  je(.  be  the  space  of  all  these  functions. 

This  can  be  iterated:  The  local  mappings  also  define  a global  mapping  from 
the  base  mesh  to  the  mesh  M'  by  HN°''  = H'~l  o • • • o HN° , for  each  i = 
Nq  + 1, . . . , N.  This  allows  to  define  the  approximation  spaces  K,'  from  the 
upper  approximation  spaces  as 

K?  = {/ o HNo'1  |/€fc*}C  ?(MNo). 

The  second  step  is  to  define  the  operators  P'  : K>+1  — > K' . Fortunately, 
working  directly  in  the  approximation  spaces  is  not  required:  Because  they  are 
isomorphic  (by  construction)  to  the  upper  approximation  spaces,  it  suffices  to 
define  some  operators  P'  : £'+1"  — > £»,  and  for  the  purpose  of  visualization, 
the  least-square  projection  operators  will  be  used.  The  operator  P'  is  thus 
defined  by 

P'  : fi+ 1 = ft+ 1 o HN°’i+ 1 Pi(/’+11’<)  o H" 

fi+i  € K,+l  defines  to  which  we  apply  one  decomposition  step  to  get 

an  approximation  /*  € K'  and  detail  coefficients  (1  in  the  linear  case,  and  2 
in  the  constant  case).  The  entire  operation  is  then  repeatedly  applied  to  /,, 
fi- 1,  etc.  Consequently,  in  practice,  everything  happens  in  the  upper  approx- 
imation spaces  which  have  a much  simpler  structure  than  the  corresponding 
approximation  spaces. 

4.3.  Results  and  remarks 

The  output  of  this  algorithm  is  a coarse  function  /yv0  defined  on  the  base 
mesh  and  a list  of  detail  coefficients  that  allow  the  exact  reconstruction  of 
the  original  function  through  the  hierarchy  of  LODs.  Figure  3 shows  some  re- 
sults in  the  linear  setting.  In  these  examples,  the  geometric  criterion  guiding 
the  decimation  priority-queue  is  simply  the  distance  from  the  1-neighbours  of 
a candidate  vertex  to  their  least-square  approximation  plane.  On  the  upper 
right  snapshot,  we  see  the  drawback  of  a geometric-only  driven  priority  queue: 
Some  quasi-planar  areas  on  the  object  have  been  severely  decimated,  leading 
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Original:  32000  vertices  Geometric  criterion:  3000  v. 


Geom.  & data  criterion:  3000  v.  Same  without  vertices 

Fig.  3.  Bracket:  linear  approximation  examples. 

to  a quite  coarse  approximation  of  the  function.  The  lower  snapshots  show 
the  result  using  the  geometric(|)  and  data(|)  based  criterion.  The  resulting 
approximations  are  better,  but  the  geometry  presents  some  visible  deforma- 
tions (this  is  also  partially  due  to  our  simple  geometric  criterion).  Finding  a 
compromise  in  an  automated  way  seems  to  be  a difficult  task.  Moreover,  if  the 
simplification  is  just  a process  prior  to  other  computations,  such  as  can  be  the 
case,  e.g.  in  mechanics,  then  a high  accuracy  in  the  approximated  function 
might  be  the  primary  interest.  Thus,  it  seems  better  to  let  the  weights  depend 
on  the  application,  under  user  control. 

4.4.  Comparison  to  classical  MRA 

In  [6],  a MRA  for  subdivision  surfaces  is  used  to  handle  both  the  geometry 
and  the  scalar  attributes  of  a mesh.  The  presentation  given  above  makes 
the  link  between  these  methods  and  the  decimation  approach.  Indeed,  from 
the  “upper”  point  of  view  — the  decomposition  using  the  operators  Pl  — it 
compares  to  decimation  in  many  respects,  whereas  it  is  also  a decimation  step 
corresponding  to  the  approximation  spaces,  which  is  exactly  what  is  done  is 
[6]  in  the  nested  case.  However,  parameterizations  are  then  obtained  without 
a local  projection  hypothesis,  thanks  to  the  particular  l-to-4  splitting  strategy 
that  is  performed  on  the  base  mesh. 
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